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SUMMARY 


This report contains the research on the CFD code evaluation with 
emphasis on supercomputing in reacting flows. Advantages of unstructured 
grids, multigrids, adaptive methods, improved flow solvers, vector processing, 
parallel processing, reduction of memory requirements are discussed. As 
examples, we include applications of supercomputing to reacting flow 
Navier-Stokes equations including shock waves and turbulence and combustion 
instability problems associated with solid and liquid propellants. This report 
does not include evaluation of codes developed by other organizations. Instead, 
the basic criteria for accuracy and efficiency have been established, and some 
applications on rocket combustion have been made. Research toward an ultimate 
goal, the most accurate and efficient CFD code, is in progress and will continue 

for years to come. 
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1. INTRODUCTION 


The literature on computational fluid dynamics is preponderant, interwoven with 
success and failure. Stable and accurate solutions are attributed to low Reynolds number, 
low Mach number, and low Pechlet number. As these numbers increase, however, we 
encounter convection of flow to be critical, resulting in turbulence for large Reynolds 
numbers, and thermal discontinuities for large Mach numbers, and thermal discontinuities 
for large Pechlet numbers. These physical phenomena lead to unstable, nonconvergent, 
and inaccurate solutions. 

It is, therefore, important to devise the most efficient schemes to obtain stable, 
convergent, and accurate solutions for difficult flow situations such as those occurring in 
propulsion systems and the space shuttle main engine in particular. Many such attempts 
have been reported in the open literature for the past two decades. No consolidated efforts, 
however, have been made to evaluate all of the available codes to date. Each code has 
merits and demerits. Often defects are buried and unknown to the user, a cause for great 
frustration. 

If the fastest time-scale in the system of equations is not physically relevant for the 
problem at hand, much lower modes carry the physical information as exhibited in implicit 
methods [1—4]. A typical example is the viscous flow past an airfoil where the speed of 
sound limitation in the cells covering the boundary layer would impose severe time step 
restrictions for an explicit scheme without achieving any further accuracy. Similarly, 
widely disparate time and length scales which occur in turbulent flows or highly convective 
flows will require special treatments. 

The use of the finite difference method (FDM), both implicit (e.g., Beam and 
Warming [3]) and explicit (e.g., McCormick [5]) has been successfully applied to 
aerodynamic problems. Initially, most applications were for external flow problems. More 
current applications have shown these methods to be useful also for internal flow problems. 
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Other approaches which are becoming useful are the finite element methods (FEM) 
using unstructured grids. Computational difficulties involved in flows with high Reynolds 
numbers, Mach numbers, and Pechlet numbers are treated in FEM [7] in a manner similar 
to FDM [8] in some cases, but there are unique features of FEM which warrant further 
investigation. These include efficient treatments of widely disparate length and time scales 
in convection, shocks, turbulence, and reacting fluids through multigrid adaptive methods 
in supercomputing with vector and parallel processing. 

For many years CFD grew around FDM as they were simple to understand and 
code, easy to vectorize, in structured grids for simple geometries. However, as computers 
became bigger and faster, attempts were made to simulate more and more complex flow 
domains, and it appeared that unstructured grids may be flexible enough to describe these 
domains. It was at this point in time that unstructured grids with finite elements - a 
natural way of discretizing operators on them — entered the scene of CFD. Since then, 
unstructured grids have been used in FEM together with domain splitting [9] and adaptive 
refinement [10, 11]. Many more developments, however, are still needed in order to 
transform FEM into efficient engineering tools for CFD. 

The Beam and Warming scheme, in its basic form, can be rewritten for 
unstructured grids. However, the approximate factorization used for structured grids must 
be replaced by the solution of a full matrix [4]. The solution of full matrices can only be 
attacked via unstructured multigrid methods. MacCormick’s implicit two-step procedure 
makes heavy use of upwind-differencing, thus always assuming a structured grid, and for 
this reason cannot be used in the present context. A barely implicit scheme [6], unlike the 
two former methods, treats only the sound waves implicitly by solving a modified Poisson 
equation for the pressure. Therefore, instead of solving five coupled equations in 3— D for 
the Euler equations, only one needs to be solved. The resulting Poisson equation is again 
solved via unstructured multigrid methods. 
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We shall discuss unstructured grids, multigrid methods, adaptive methods, 
improved flow solvers, supercomputer utilization, applications to turbulent reacting flows, 
solutions using finite element techniques, and applications toward combustion instability in 
the following sections. 

2. UNSTRUCTURED GRIDS 

The accurate representation of arbitrary domains represents perhaps one of the most 
challenging problems in CFD. The magnitude of this problems does not become apparent 
in two dimensions because only a few singular points usually appear in the field (and may 
be ignored) and, due to the computer capacity now available, a gross overmeshing in 
certain regions of the domain can still be handled. 

However, anyone trying to mesh complicated geometries in three dimensions with 
structured grids will encounter singular lines, and the unavoidable cost of overmeshing can 
no longer be ignored (the result being coarse grids). 

It is by now generally accepted that unstructured grids are capable of describing 
accurately complicated geometries in 3— D better than structured grids. Two different 
levels of unstructuring are possible. 

(1) Macro-unstructuring, where blocks of structured grids are combined to form an 
unstructured grid (these are so-called zonal methods). 

(2) Micro— unstructuring, in which case the point and element distribution can, in 
principle, be random. 

Although macro-unstructuring is being actively pursued by several groups [12—14], 
the inherent structure at the difference— level precludes simple mesh refinement by local 
enrichment (which destroys the grid structure). Major problems will also appear when the 
region to be refined/enriched crosses zone— boundaries. Micro— unstructuring does not have 
these inherent limitations, but does have disadvantages such as programming difficulties 
and large storage requirements. 
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However great the disadvantages of micro— unstructured grids may appear, the 
advantages these grids offer by far outweigh them. For example, (1) any geometry can be 
described [15,16], (2) mesh refinement either by movement [17—20], enrichment [21—31], or 
remeshing [32] presents no problems, and (3) domain splitting [9, 46] for transient problems 
can easily be performed. Recently, Jameson [16] also presented results using unstructured 
grids, and this alone may indicate a turning point in the development of CFD. 

The fast generation of grids for arbitrary domains in three dimensions has been the 
focus of much research in recent years. A variety of different approaches have been 
investigated. The most promising seems to be: the macro— element approach [34, 35], 

Watson’s algorithm [36 — 41] for Voronoi tesselations combined with a point distribution 
obtained by superposition of local (structured) grids [16], modified octree [42], on an 
advancing front [32,43,44], and from a regular background [45]. All of them have 
advantages and disadvantages, and none has been fast and simple enough to generate 
efficiently grids of the size needed in 3— D aerodynamic simulations. Grid generation for 
unstructured grids in 3— D has been pursued for only a few years by a few individuals. 

Based on the current state— of— the— art, the following conclusions are drawn: 

• Partially unstructured grids (the macro— element approach [34,35]) do not offer 
enough flexibility to serve as the basis of a general mesh generator for complex 
domains, unless some major breakthrough in interactive graphical display is 
achieved. 

• Use of regular background grids [45] is not advisable for unbounded problems, as 
element stretching and point clustering are not accommodated easily. However, this 
technique may prove useful for internal flow problems. 

• The generation of points via superposition of local grids (mapping) [16] seems to 
offer the greatest flexibility at minimal cost. The point distribution for each local 
grid is obtained algebraically and is therefore very fast. The grids are chosen from a 
menu of possible local grids. 
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• The tetrahedrization of the domain via Watson’s algorithm [36,37], as employed and 
modified in [16,18,41], is not advisable as this algorithm is suboptimal, requiring 
0(N 1 ' 5 ) operations, and the treatment of voids in the fluid domain becomes both 
grid logic and CPU— intensive. 

• It seems attractive to pursue the element generation of [44], combining it with a fast 
neighbor finder [46—48]. 

3. MULTIGRID METHODS 

Multigrid methods combine two very desirable properties in that they require the 
least amount of operations to solve large problems (0 (N logN) for a problems with N grid 
points) and their storage requirements are also low (again, 0(N logN) for a problem with N 
grid points). In 1985, Lohner and Morgan [46, 49] advanced the concept of unstructured 
multigrid methods. It became clear that as the finest grid had to be unstructured in order 
to accurately represent the domain, it could not be obtained by subdivision of some coarser 
grid. Instead, a set of unrelated coarsening grids had to be employed. The reason why 
multigrid methods should still work on sets of unrelated, unstructured grids — the same 
argument on which all multigrid methods base the convergence rate estimates — is that if 
the residual is smooth, any coarser grid should be able to "see" it. In all other aspects, the 
theory follows exactly the lines of traditional multigrid methods [50,51]. 

The solution of elliptic PDE’s via multigrid methods is by now well understood, and 
rigorous theoretical estimates for the expected convergence rates are available [51]. The 
main difficulty that can appear for unstructured grid lies in the construction of efficient 
smoothers, as neither line— nor plane— relaxation are possible. If Jacobi— type smoothers 
are employed, the convergence rate of the highest modes can degrade seriously for highly 
stretched elements or diffusion tensors in which one direction is dominant [52]. Three 
different smoothing schemes are known to avoid this problem: 
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(1) Use of supersteps: here, the simple Jacobi— smoothing is over— and under— relaxed 
alternately using a Chebyshev— series. Although not advisable for highly stretched 
grids, e.g., stretching beyond 1:100, this method is very simple to code and lends 
itself easily to vectorization. 

(2) Solution of local problems: instead of a tighter coupling of modes via line— or 
plane-relaxation, groups of elements are relaxed, producing the desired effect [53]. 
This method is applicable in all cases, but may not be vectorizable and also requires 
some software— overhead. 

(3) Element-by-element preconditioning: although the transfer of information in the 
element— by— element iterative solver [54,55] is local in nature and therefore cannot 
compete with multigrid methods, this scheme may prove useful as a preconditioner. 
The compression of the eigenvalue spectrum is achieved by multiplying the system 
matrix with the (local) element— matrix inverses where appropriate. Vectorization 
of this type of method should also be investigated further. 

For the hyperbolic case, the theory of multigrid methods is still far from complete. 
Although Ni’s method [56—59] has been shown to work well in many cases, Jameson’s 
multigrid solver [60,61] seems to emerge as the more reliable. This is to be expected, since 
the Runge— Kutta time stepping allows more possibilities for choosing appropriate 
"damping— sequences" than the Lax— Wendroff schemes. The combination of unstructured 
multigrid methods with Runge— Kutta time stepping for the Euler equations also appears to 
be useful [62,63]. 

4. ADAPTIVE METHODS 

In convection— dominated problems, discontinuities or regions with sharp gradients 
will usually appear. The regions in which the flow variables vary abruptly are usually 
small and are surrounded by large portions of the field in which the flow varies smoothly. 
Therefore, it is attractive to locally and adaptively refine the mesh where needed until a 
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preset tolerance for the error has been achieved. Because of the obvious advantages of 
adaptive refinement, this field is currently receiving increased attention in the literature 
[10,11]. Any adaptive refinement scheme consists of three different stages: determining 
what we want to achieve by refining the grid, developing an error indicator/estimator to 
detect the regions to be refined, and a refinement strategy such as movement, enrichment, 
or remeshing. There are three specific directions to be considered: 

(a) Typically, one aims to have an equidistribution of the "error" throughout the grid 
[10,11,51-63]. 

(b) Whole families of error indicators based on different concepts have been shown to be 
useful. Among the most popular are those based on the change of some "indicator 
variable" (e.g., entropy [27] or Mach number [20]), those based on interpolation 
theory estimates [18, 21-26,64,65], and the indicators based on Richardson 
extrapolation [66]. 

(c) For the accurate resolution of the flowfields at or near discontinuities, 
P-enrichment does not seem to be attractive (besides, P-enrichment implies a 
considerable increase in software complexity). However, it may prove useful for 
boundary layers, where an essentially smooth flowfield needs to be resolved. 

Although further research on more reliable error indicators [65] is needed, the 

directional refinement advanced in [31] appears to be satisfactory. It is based on the 
observation that in most flowfields the regions that ought to be refined are of lower 
dimensionality than the physical space in which the solution is sought. Therefore, if thin, 
elongated elements parallel to those discontinuities could be generated during the adaptive 
refinement process, considerable savings in CPU and storage would be realized without 
sacrificing accuracy. 

The first algorithm devised for this purpose was based on mesh enrichment [31], and 
turned out to be storage, CPU, and software intensive. Since then, Palmerio and Dervieux 
[20] have tried to incorporate this concept into a mesh movement framework, while Peraire 
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et al, [32], have advanced the concept of refinement by remeshing. This last concept 
represents a new and powerful refinement strategy, combining in a very elegant way the 
advantages of mesh enrichment (such as versatility by the introduction of points and a 
coarse initial grid for steady state problems) and mesh movement (which produces the 
desired element shapes near shocks). If it proves useful in 3— D (the search— problem needs 
to be addressed here), it may completely replace movement and enrichment as refinement 
strategies. Another major area of applications for fast remeshing algorithms is given by 
3— D Free— Lagrange Methods [66], where restructuring of the grid currently represents a 
major problem. Remeshing would be an ideal solution in this case. 

When solving transient problems in which only a few discontinuities appear, 
adaptive refinement can also be useful in reducing storage and CPU requirements. 
However, in comparison to steady state problems, further constraints need to be placed on 
the refinement algorithms in order to realize significant savings: 

(a) As the grid adaptation has to be performed many times, the adaptation algorithm 
must be fast, and therefore must lend itself to vectorization/parallelization. 

(b) As the grid adaptation process becomes an integral part of any code, the algorithm 
should not be storage intensive. 

(c) As the feature that has been refined may pass again (e.g., a shock reflection), the 
original grid should be recovered after the feature has passed. 

Of course, directional remeshing or even movement could be incorporated for those 
cases in which the discontinuities are fairly straight (for curved shocks that interact with 
each other, only classic h— enrichment will yield acceptable solutions), but the interpolation 
problem must be solved (otherwise the shock— speeds will be wrong). Another potential 
problem may arise due to the (apparent) non-vectorizability of the remeshing algorithm. 
Obviously, this whole topic of adaptive refinement for transient problems represents one of 
the most difficult ones in CFD; many further innovations are expected. 
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5. IMPROVED FLOW SOLVERS 

For unstructured grids, the extension of schemes from 1-D to 2-D/3-D cannot be 
performed by operator splitting. This means that only schemes that are truly multi- 
dimensional in nature can be used. Toward this end, Petrov— Galerkin methods [67], the 
flux corrected transport (FCT) algorithms, and those schemes generalized to multi- 
dimensional problems [71-74] appear to be most suitable. 

The idea behind FCT is to combine a high order scheme with a low order scheme in 
such a way that in regions where the variables under consideration vary smoothly, the high 
order scheme is employed, whereas the low order scheme is favored in those regions where 
the variables vary abruptly. In this scheme, no nonphysical over/undershoots are created. 

It is at this point that a further constraint, given by the conservation equations 
themselves, must be taken into account: strict conservation on the discrete level should be 
maintained. The simplest way to guarantee this for the node-centered schemes considered 
here is by constructing schemes for which the sum of the contributions of each individual 
element (cell) to its surrounding nodes vanishes, i.e., "all that comes in goes out". This 
means that the limiting will be carried out in elements (cells). 

There are several aspects in FCT worthy of mention: 

(a) For systems of equations, no obvious or natural way of limiting has been identified 
yet. Several possibilities have been explored, among them (for the Euler equations) 
operator splitting (treating each equation independently), the use of the average of 
the limitors for each equation, limiting based on some "key variable" (pressure, 
entropy), and others [74], 

(b) Entropy is not always monotonic for FCT. This may be due to the low order 
scheme employed. It is obvious that the ultimate low order scheme is Godunov’s 
scheme [75—78], but this scheme is much more expensive that the simple smoothers 
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currently in use [73,74]. A simpler, more efficient version of Godunov’s scheme for 
unstructured grids should be developed. 

(c) Steepeners for contact discontinuities: as contact discontinuities are linear, any 
scheme that does not possess a steepener will eventually flatten these out. Note 
that no physical argument leads to a distinction of linear and nonlinear 
discontinuities. Some kind of detection mechanism may be usefully incorporated for 
contact discontinuities. 

6. SUPERCOMPUTER UTILIZATION 

However good a method may be, if it does not lend itself to some form of 
parallelism, its future will always remain a dubious one. The speed-up ratio between a 
code that exploits the machine hardware and one that does not lies between 1:10 and 1:20 
on today’s vector machines. This performance ratio will go up drastically when massively 
parallel computing becomes available [79,80]. Fortunately, the bigger the problem to be 
solved, the easier it is to exploit some inherent parallelism of an algorithm. 

6.1 VECTOR PROCESSING 

On vector machines the important factors that determine the performance of an 
algorithm are DO-loop length and contiguity in memory (even on a CRAY!). Normally, 
the vector length in CFD codes for element subroutines is of order 8—10, far too short to 
exploit vector— machine hardware. The only way to achieve acceptable vector— lengths is to 
perform the assembly process on groups of elements, possibly the whole grid at once. This 
means that one— element type codes should be favored, in contrast to the more usual 
many— element type code now in use in industry. 

Three different types of DO— loops are most often encountered: 

(a) Loops over the same type of data: these are loops which involve only one type of 
data (either point or element data) and are the "favorites" of vector machines. 
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(b) GATHER— loops: this type of loop appears when point information needs to be 
processed at the element level (a point may be shared by several elements). Most 
vector machines have hardware-tailored GATHER devices, but this type of loop 
will, nevertheless, run between 2.5-5. 0 times slower than type (a). 

(c) SCATTER— ADD: these loops occur when assembling element contributions at 

points (e.g., formation of a right-hand side vector). As a point may receive 
contributions from several elements in the same loop, the simple SCATTER 
operation is not sufficient [81—83]. Coloring schemes have to be devised, and the 
original loop (over all the elements) has to be broken up into groups of elements so 
that the use of straight SCATTER operations becomes possible. This type of loop 
will run at roughly the same speed as the GATHER-loop. 

6.2 PARALLEL PROCESSING 

When speculating about parallel computing, one ought to distinguish between 
mildly parallel machines (up to 10 processors) and massively parallel machines. Examples 
of mildly parallel machines are the CRAY XMP/48 and the systems of array processors 
attached to a host machine. From a user’s point of view, the main difficulties facing the 
development of codes on such a system are the poor FORTRAN capabilities and the very 
small local memory of the presently available array processors, as well as the very bad 
debugging software. Therefore, it pays off the rewrite only codes which have been 
thoroughly tested and will undergo no major modifications for these machines. In an area 
as dynamic as CFD, few codes ever made it to an array processor. However, FCT models 
have and with great success [84,85]. 

From a theoretical point of view, all algorithms which split the domain into 
sufficiently large subdomains are suited for mildly parallel machines. Examples of this 
kind are simple operator splitting (e.g., line by line, as long as the line is big enough), or 
the growing family of domain splitting algorithms [9,33]. 



12 


6.3 REDUCTION OF MEMORY REQUIREMENTS 

Schemes operating efficiently on unstructured grids require more memory than their 
structured grid counterparts. Among the possibilities that can be pursued in order to 
reduce memory requirements, we mention the following: 

(a) Careful coding: an obvious possibility, but one that is usually considered only after a 
code has been shown to work (i.e., at a stage where as little as possible should be 
changed). Experience indicates that 30% reduction of memory may be achieved at 
the expense of 20% increase in CPU. 

(b) Splitting into subdomains: the peak efficiency of vector machines is achieved for 
vector lengths that are only a fraction of the total number of grid points typically 
encountered when solving 3-D problems. The idea is then to save as much storage 
on temporal arrays as possible, performing all algorithmic steps (formation of 
right-hand sides, limiting, etc.) on subdomains. Memory requirements may be 
reduced by about 50% at the expense of some additional CPU time. 

(c) Encasement of the unstructured grid in a structured grid: here, unstructuredness is 
only allowed close to the body where the highest geometrical/physical complexity is 
expected. At wider distances from the body, a structured grid is employed. This 
approach holds considerable promise of reducing requirements by more than an 
order of magnitude, but will require sophisticated programing and mesh generation 
capabilities. 

7. SUPERCOMPUTER APPLICATIONS TO ROCKET COMBUSTION 
INSTABILITIES 

7.1 GENERAL 

The growth or decay of energy is responsible for instability or stability of waves in 
fluids, respectively. In general, there are three different sources of energy growth or decay. 
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They are: (1) pressure fluctuations (acoustic waves); (2) velocity fluctuations 

(hydrodynamic waves due to transition to turbulence, shear layers, or shedding of 
vortices);and (3) density fluctuations (intrinsic waves due to compressibility and/or chain 
reactions of unstable chemical radicals in combustion). Extensive research works have 
been carried out for acoustic instability [135-140] and hydrodynamic instability [140-143], 
with limited studies available on intrinsic instability [144]. Flandro [145] presented the 
energy balance method in which the acoustic energy equation was used to develop the 
nonlinear stability equation based on isentropic assumption and linear superposition of 
entropy change. Each of the above types of instability requires different methods of 
analyses and unification of such methods has not been attempted. The objective of the 
present paper is directed toward our desire to introduce a unified general method for wave 
instability analyses. 

The basic approach toward this achievement, called the entropy controlled 
instability (ECI) method, is founded on the concept of entropy changes in which 
nonlinearity is a prime factor. This is because our ultimate goal is a general theory which 
enables nonlinear waves to be studied in detail without making simplified initial 
assumptions such as isentropy. Nonlinear waves invoke nonisentropy in which growth or 
decay of energy is properly accounted for. This line of thinking is dictated also by our 
desire to deal with compressible fluids, with a special case reduced to incompressible flows. 
To this end, we examine the energy equation written in terms of entropy changes as well as 
fluctuations of pressure, velocity, and density. Introducing spatial averages, or applying 
the Green— Gauss theorem, we integrate the energy equation by parts through which 
boundary surface integrals arise, playing the role of acoustic intensities. Additional terms 
representing the spatial growth of energy also arise in this process of integration by parts. 
Performing the time averages, we finally arrive at the nonlinear, nonisentropic stability 
equation, which assumes the form of nonlinear integro-ordinary differential equation for 
the energy growth factor indicative of stability or instability of waves. The integrands of 
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this nonlinear stability equation (the word ‘nonisentropic’ omitted but implied hereafter) 
are contributed from the results of the Navier-Stokes solutions using the Taylor-Galerkin 
finite elements [146—150]. The nonlinear stability equation will then be solved using the 
fourth order Runge— Kutta method to calculate the energy growth factor. It is shown that 
the general theory presented here is reduced to a simple linear theory of Cantrell and Hart 
[135] or Culick [136] as a special case with appropriate assumptions. 

We shall discuss the details of these processes in the following sections. 

7.2. ENTROPY CONTROLLED INSTABILITY (ECI) METHOD 
7.2.1 NAVIER-STOKES SOLUTIONS 


In order to perform stability analysis, it is necessary to first solve the unsteady 
Navier-Stokes equations representing the desired physics. To this end, the most general 
conservation form of Navier Stokes equations for compressible flows is written as 
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where is the viscous stress tensor 
r ij — P ( v i>j + V j>i — J v k>k^ij) 


( 1 ) 


(la,b) 


(lc,d) 


(le) 



15 


and E is the stagnation energy 

E = e + \ v i v i = c p T “ p + k v i v i (If) 

and f ki is the body force and qj is the heat flux vector. 

qj — — AT,j + p D ^ H k Y k ,j (lg) 

Here, A and D are the thermal conductivity and mass diffusivity, respectively. H k is the 
total enthalpy of species k, Y k is the mass fraction for the species k, and w k is the reaction 
rate for the species k. Note that all derivatives are covariant in case of cylindrical 
coordinates. 

For turbulent flows we add to (1) additional transport equations for turbulent 
kinetic energy and dissipation energy (k— e model), respectively, 

It W + ^ w) = G ~ pe ( lh ) 

VpM + liE^V ~ = c i £ G ~ c 2^ e V k (li) 

where G is the turbulent thermal dissipation energy transport and 

M t Jr 2 

M k = M + ^ , M = P%~ e , 

c^ = 0.09, q = 1.44, c 2 = 1.92, a k = 1.0, <r e = 1.3 

Appropriate modifications to the equations of momentum and energy and the suitable 
turbulent combustion model should be included as detailed in [146]. 

The solution procedure using the Taylor-Galerkin method are found in [147-151], 
which will not be repeated here. Examples include shock waves contributing to nonlinear, 
nonisentropic combustion instabilities for the space shuttle main engine combustion/thrust 
chamber and the side— burning solid propellant rocket motors. Calculations of energy 
growth factors and stability criteria are discussed below. 
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7.2.2 NONLINEAR, NONISENTROPIC STABILITY ANALYSIS 


With the solutions of unsteady Navier— Stokes equations for the time— dependent 
periodic oscillatory initial boundary conditions available, we now turn to the entropy 
controlled instability (ECI) method of determining the unstable wave phenomena. To 
introduce entropy changes in the energy equation, we begin with thermodynamic relations 
for an ideal nonisentropic gas and write the pressure gradients in terms of density and 
entropy gradients, 


P,i - fl2 Pn + %r- S,i 


( 2 ) 


where a is the speed of sound, S is the specific entropy, and the comma denotes the partial 
derivative with respect to Xj. The gradient of the stagnation energy, E, assumes the form 

E, i = (c p T — p + j v j v j)>i (3) 

where the repeated indices imply summing. Performing the differentiation implied in 
(2) results in 


C V C v p 

E ’i = K^ P >i T j ^>i + Vj’i 

P 

where R is the specific gas constant. Substituting (2) into (4) yields 
P E >i = £p.i + ft S.i + pVjVj.j 

Let us now examine the conservation form of the energy equation 

d 


3t 


(pE) + (pE Vi - (Ty Vj + qj.i^i = 0 


where a y is the stress tensor 

^ = - ptfy + n (v ifj + Vj,i - 1 v k)k 
Substituting (5) into (6) gives 


(4) 


(5) 


( 6 ) 


(7) 


It (^ E ) + v i)>i + v i[p Pn + ft S >i + P v j v j.i] - (^ijVj),i + Qi.i = 0 (8) 

This is the entropy controlled energy equation, instrumental in determining the nonlinear 
wave instability. 
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The pressure, velocity, and density may be written in the form 


p = p + p' 

(9a) 

Vi = Vj + v[ 

(9b) 

p = p + p' 

(9c) 


where the symbols, bar and prime, denote the mean and fluctuation parts, respectively. 
From thermodynamic relations we may write the entropy difference in the form 



Expanding the RHS of (10) in infinite series and after some algebra, we obtain 

S = R[S(1) + S( 2 ) + S( 3 ) + S( 4 ) + ...j + S 0 (11) 

where S 0 is the entropy at the initial state (S 0 = 0), and 



Here the terms with fifth order or higher are neglected, assuming that they are negligible. 

To obtain acoustic energy equation, we apply the Green— Gauss theorem or integrate 
(8) by parts and take a time average in the form 


<} n dQ + J r E ^ v i n i + v i(P n i + £ S n i + P v j v j n i) “ *ij v j n i “ A ( T ) T »i n i 
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+ P D k ?i C pk T Y k'i “■ dF ~ } n [ E >i P V i + f ( v i ^).i + ( v i f).i S + Vj 

+ (^T)),! T„ — (pD k E c pk T)„ Y klI ]dn ) = 0 ( 12 ) 

where Q and T represent the domain and boundary surface, respectively, ( • ) implies time 
averages and n i denotes the component of a vector normal to the surface. 

It is interesting to note that the domain integral terms (last four terms) arise as a 
result of the integration by parts implying the spatial growth of energy. This spatial 
growth of energy is in addition to the usual temporal growth of energy indicated by the 
domain integral with the time derivative (first term in (8)). Let e be the energy growth 
factor, stable for 0 < e < 1, unstable for e > 1, with e=l indicating the neutral stability. 
The fluctuation terms in (9) and (11) contribute to (12) as multiples of higher order 
fluctuations. Multiplying the like powers of e with the fluctuation terms of the same order 
and retaining the terms up to and including the fourth order, we obtain 



[ fl (*0 C ^1 + f2 ^2 + f3 ^3 + f4 ^4 + •■•) dfl) 

(13a) 

| r (-) dr > = < 

+ f3 ^3 + £ V 4 + •••) dr) 

(13b) 


Here 6 0 and <p 0 contain only the mean quantity, 6 V 6 2 , 6 V S 4 and <f> v <p 2 , <p v <p 4 are referred 
to as the first, second, third, and fourth order perturbations capable of growth and decay in 
acoustic energy as dictated by the magnitudes of e, growing if e > 1 and decaying if e < 1. 
Notice that the relations in (9a, b, and c) imply that the Navier— Stokes solutions are 
obtained with the initial condition 6=1 and the energy in (13) may grow or decay in 
accordance with the magnitude of e from the initial or reference state 6 = 1. 

Performing the differentiation as implied in (12), we obtain 

(f 2 E l + 6 3 E 2 + 6 4 E 3 ) = 6 2 I t + 6 3 I 2 + 6 4 I 3 (14) 

in which the zeroth order terms are canceled and first order terms vanish due to time 
averages. Notice that the time averages are contained in E ( i} and I ( f) . Here, the energy 
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growth factor e is an explicit function of time. However, E (i) is no longer an explicit 
function of time because of its time averages. Thus the partial derivative with respect to 
time in (14) involves only e, not E( ^ , so that 

d c f2I i + f3l 2 + e<1 3 

^ 2£E t + 3 c 2 E 2 + 4r 3 E 3 

or 

A if ^^2 TO 2 Esi ^ 

^ = (£ii+c 2i 2 +£ 3i 3 ) 2 ^{ i -e 2E ^+ (f;) _ et"]} ( 15 ) 

where higher order terms and those terms much smaller than unity have been neglected. 

It follows from (15) that the nonlinear stability equation takes the form 

^ — aj£ — a 2 e 2 — a 3 € 3 = 0 (16) 

This is a form identical to Flandro [145], although the basic approach to the formulation 
and the solution procedures differ. Here a t , a 2 , and a 3 are energy growth rate parameters 
of first, second, and third order, respectively. 




(17a) 

1 3E 2 
a 2 = set; ( ! 2 - ! i) 


(17b) 

1 It 3E2 t 4- \ 9 (\ 2E ’l 1 1 

a 3 = 5T^{ I 3-^ 1 2+ E7J hj 

(17c) 

with 



E -<J/ 1 ’ dn > 


(18a) 

E,= ([ a«» dfl) 


(18b) 

E, = <j n a'» ) dfl) 


(18c) 

I, = <[ 6 (l1 dfl) — ( 
J Q * 

c ( ii> mdr) 

r 

(19a) 
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I, = <j n 6< !’ dll) - <J r c)” n, dr) 

(19b) 

I, = <J n 6' »> dtl) - <j r ci«> Hi dr> 

(19c) 


Explicit forms of a ( 11 , i» ( 1} , and c ( l) will be discussed in Section 3 and the integrands for 
the second and third order growth parameters are shown in Appendix A. The physical 
significance and interpretations of the above process will be discussed in the following 
subsection. 

7.2.3 PHYSICAL INTERPRETATION 

The physical significance in the above development may be schematically 
demonstrated as shown in Fig. 1. First, unsteady Navier— Stokes solutions with time 
dependent oscillatory initial boundary conditions are obtained at each computational grid 
as shown in Fig. la. A typical variable vs time, (Fig. lb), at any grid point (finite element 
node) may exhibit sawtooth type wave forms as well as sinusoidal waves. Imagine that 
several variables at each of the several thousand nodes present themselves in the form as 
shown in Fig. lb. Each wave form of different frequencies at different spatial locations is 
expected to contribute to the overall stability or instability. Take one complete peak wave 
form period as indicated by r = nAt where At is the Navier— Stokes time step with n 
usually being in the range between teens and hundreds. 

Performing the time average for this time period (nAt), we obtain the mean 
quantity (f) of the variable f. Then the disturbance (fluctuation) part f' is calculated as 
{' = f-f 

where f is the Navier— Stokes solution. 

From these data, spatial integrations involving both domain and boundary surfaces 

and time averages (Fig. la) as dictated by (18) and (19) are then carried out. Here, 

a ( n) dQ denotes the temporal growth of energy, whereas 6 ( n) dQ represents the spatial 

Jn Jn 
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growth of energy. The boundary integral j c ( . n) n £ dT indicates the acoustic intensity. 

For combustion chamber applications, the burning surface admittance, response functions, 
or information on nozzle entrance or acoustic liner, etc., can be initially imposed on the 
Navier-Stokes solutions as boundary conditions to generate the oscillatory motions as 
depicted in Fig. lb. These boundary data, however, reappear as called for in the integrand, 

c( i n) , * n (I®)* They act as acoustic intensity and eventually contribute to the 
determination of stability or instability of the system. 

The ingredients of a (n) , b (n) , and c ( . n) arise from fluctuation parts of the 

Navier-Stokes solutions as well as the mean parts. These fluctuations through the 

temporal growth of energy a (n) , spatial growth of energy 6 ( n) ) and acoustic intensity c ( n) , 

i 

are responsible for driving the system toward instability. The consequences lead to 
determination of the energy growth factor vs time as shown in Fig. lc by solving the 
nonlinear stability equation (16). It is important to realize that the Navier-Stokes 
solutions and the nonlinear stability equation (16) encompass acoustic waves, 
hydrodynamic (vortical or shear layer) waves, intrinsic (density fluctuations or chemically 
reacting combustion) waves, or combined effects of all types of wave interactions. 

Notice that, in (6) and (12), heat flux changes are not directly involved. This is 
because the spatial integration of fluctuations of heat flux or temperature vanishes upon 
time averages. However, fluctuations of temperature and chemical species mass fractions 
have been included in the Navier-Stokes solutions of (1), thus indirectly contributing to 
the fluctuations of density (/?') in the stability calculations. 

Some comments on the mean pressure changes (DC shifts), pressure coupling, 
velocity coupling, limit cycles, and triggering are in order. These phenomena often 
observed in solid propellant combustion chambers, would appear in the flowfield under 
appropriate initial and boundary conditions during the Navier-Stokes calculations. The 
purpose of the ECI method is merely to determine whether the flowfield containing such 
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physical phenomena is stable or unstable. It is emphasized that all physical phenomena are 
expected to prevail and be identified in the solution of Navier— Stokes equations, which are 
then reconfirmed in the stability analysis. 

7.3. LINEAR INSTABILITY 

To gain an insight into the solution of (16), we may neglect the last two terms 
associated with the second and third order energy growth rate parameters and write 

^-a t e = 0 (20) 

which yields a solution in the form 

4i c = a t t + c t (21) 

To establish an initial condition, we assume neutral stability e = 1 at t = 0. This gives c t 
= 0. Thus, the solution takes the form 

OL-t 

£ = e 1 (22) 

Under this initial condition, there exists a unique solution for any given with t > 0. The 
stability criteria are: 

Stable: 0 < e < 1 with — m < a y < 0 

Neutral stability: e = 1 with aj = 0 

Unstable: 1 < e < ® with 0 < a { < a> 

Although these criteria are not applicable for the nonlinear equation (16) when a 2 and a 3 
are involved, the same initial condition, = 1 at t = 0, as originally defined in (13a, b), 
can be used. That is, there exists a unique solution e for any given a v a 2 and a 3 with t > 
0. Therefore, the criteria for stability in terms of e remain the same with various 
combinations of a v a 2 , and a r 
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7.3.1 NONISENTRIPIC CASE 


Notice that the integration involved in the growth rate parameters a 2 , and a 3 
and the subsequent solution of the nonlinear ordinary differential equation (16) are 
formidable. However, it would be informative to examine the case of linear instability 
(a 2 = 0, a 3 = 0) given by (20) but with nonisentropy. Therefore, the first order energy 
growth rate parameter takes the form, from (20) and (17a), 


_ 1 de _ 1 1 

a i-73t 


<[ b { ‘>dn> - ([ C<1> nzdr) 

Jn Jn i 


(2 a ( i> 


dfl) 


(23) 


Assuming that the fluid is inviscid (p, = 0) the explicit forms of integrals in (23) become 


j„ a '” dn = |„ ( i v i T i + ''Vi) dn 


(24a) 


[ 6 ( O dfl = [ \pw i ( ^ — I ^ h h v]vp,j + (p\[ + Vj /?’)(— ^ — 

Jn Jn 1 (^ 1)^2 (r _ 1)/S 3 ^ 1 \r-l)~P 


— ^ 


L_2 1- y.y!) . -|- yV( 2 1- ^ y.y.) ■ + 0 ( '—f— -| _ Q ’ 2 

/ ,s - 2 J V’I r T ? V jV’l ^ P \ _o ^ _ 3 ^ 

(t-1 )p (7-1 )p P P l 


* f t l » 




P Vi p Vi p V iP pv^ PV: 

+ — )>i + p( -+— ),! + (“* 


IP 


P P 2 


(t-1)p (rl)p 


)(pvi 


+ pVi) '- + w + W )(P5i),i + (?v;vi + 


+ P’Vp’i + vj (pV;Vj + f vJVj + p'VjVj),, 


(24b) 
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Jr cU> n id r = J r { ? v i( - ^ + Jgj + \ *j*j) + if v; + V)(^ 


(r-i)p 


,'2 


— ^ + VjVj) + /?’vj (—2 + l VjVj) + p’vj + ( ^ 

(7-1 V 1 1 Vi )p 2,1 HrW 


+ ~f L —) PVi + (r* 


IP 


(r-i)p (t-i)p 


-)(p y J + p'vj + Vi (pvjvj 


+ VVjVj) + vj (2^ViVj + p’ViVi) + vj pj ni dr (24c) 

Notice that the first term on RHS of (24a) represents the kinetic energy of the 

sound wave. If isentropic assumption is made, then p' can be expanded into infinite series 

in terms of p’, p’ 2 , ... such that the second term on the RHS of (24a) contains the potential 

energy of the sound wave usually identified in the linear acoustic energy equation. 

However, our objective here is to allow entropy to change and, therefore, we must keep p' 

to remain subjected to nonisentropy. Another important observation is the domain 

integral for b l l) which has appeared for the first time as a result of the integration by parts 

(24b), signifying the spatial growth of energy. The boundary integral for c ( 1] n { contains 

i 

the terms of acoustic intensity normally identified in the linear acoustics, but significantly 
in a different form. Additional terms which arise in the process developed in this analysis 
will allow determination of stability or instability of wave motions with explicit changes of 
entropy contributing to the growth or decay of energy. However, it is not possible to 
evaluate the integrals (24a, b, c) because the density fluctuations p' cannot analytically be 
determined. This difficulty can be resolved in a certain special case if isentropy is assumed 
as discussed below. 


7.3.2. ISENTROPIC CASE 

If the flow is isentropic then it can be shown that 
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where a denotes the speed of sound without flow. As far as the linear stability is concerned 
only the first two terms on the RHS of (25) will contribute to (24a, b, c) as seen from 
substitution of (25) into (9c) and subsequently to (24a, b,c). 


r r t>' 2 V:V:p 

a ( i ) dQ= (fvjvj+E — + JJ—)dn 

Jn Jn 2 J J 2pa 2 a 2 


dn 


= L {*■ ( 


1 h \ VjVj + ...), i + ...j dn 

\nW ( -v<— 1 ) 75®n ^ 2 J J J 


(?-l)p 2 a 2 {j-l)ph A 


r r f P' 2 Vi pviVjV. > 

Jr c i nidr = Jr H + + " + IT- + ~j ^ 


_ , , P^i^j v i 


(26a) 


(26b) 


(26c) 


Since p' does not appear in the integrals of (26a, b, c) it is now possible to substitute 
analytical forms of p' and vj in terms of time dependent acoustic eigenfunctions and 
perform explicit analytical integrations. However, it is dear that the additional domain 
integral (26b) and many more additional terms arising as a result of the entropy controlled 
energy equation given by (8) or (12) will produce the results quite contrary to the 
traditional isentropic solution, (see Cantrell and Hart [135] or Culick [136]). 


r r P ' 2 v i 

- < IrH + 

2 ( J n <5? v i v i + 


P v i v j v n 

+ P v iVj +— 7 — n i dr > 


f-; + ^ 

2 pa 2 a 2 


Note that the terms on the numerator are the same as the first four terms of (26c) with the 
negative sign as indicated in (19a), and the terms of the denominator are identical to those 
in (26a). However, the domain integral of (26b) is absent. This integral contributed by 
b ( represents the spatial growth or decay of energy as balanced by the boundary 
conditions and acoustic intensities on the boundary surface. This is in contrast to the 
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temporal growth or decay of energy originating from (26a). 

Despite the special feature in the proposed formulation, however, the analytical 
forms for p' and v- as used by Cantrell and Hart [1], are incapable of simulating sawtooth 
type shock waves. For example consider the acoustic field of a cylinder with the radius R 
and length L, 

p' = me [p cos ^ J m (/)„,, £) cos (mO) ei-‘] (28a) 

v; = We [ COS J„ (/3„ n i) cos (m «) e“i] (28b) 

1 UJP J 

Here angular frequency, u, is defined as 



with i, m, n, being the positive integers or zero, J^, the Bessel function of order m, and (3^ 
the nth root of 3^(0) = 0. Retaining only those terms arising from the standard acoustic 
energy equations (without employing the entropy controlled energy equation), a simple 
solution can be obtained using (28a, b) to calculate the linear energy growth rate 
parameter as demonstrated by Cantrell and Hart [135]. If the entire terms implied in 
(26a, b, c) are used, however, it is no longer possible to obtain analytical solutions. 


7.4 NONLINEAR INSTABILITY 


Nonlinear, nonisentropic waves occur in many industrial propulsion combustion 
systems. Shock waves may interact with turbulent vortical waves or shear layers of liquid 
jets in gas medium. Density fluctuations due to chemical reactions may also be combined. 
We have seen that analytical solutions of even the linear instability by ECI method 
presented in the previous section are intractable. First of all, the fluctuation variables p', 
v-, and p' are to be numerically calculated by solving the Navier— Stokes equations. Then 
we must perform numerical integrations required to evaluate the energy growth rate 
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parameters a 2 , and a 3 . Subsequently, numerical solutions of the nonlinear ordinary 
differential equation must be carried out. 

Numerical integrations as required by (18a, b, c) and (19a, b, c) can be performed 
most efficiently by Galerkin finite element techniques to calculate the energy growth rate 
parameters o t , a 2 , and a 3 according to (17a, b, c). Finally the nonlinear ordinally 
differential equation (16) is solved using the Newton— Raphson method to determine the 
energy growth factor e. Thus, it is seen that the determination of stability (0 < c < 1), 
instability (e > 1), or neutral stability (e = 1) is made available for each nAt period and 
we move on until desired time is reached as shown in Fig. lc. It is interesting to see that, 
in this decision making process of stability or instability, "all" nodal points (thousands of 
nodes) with their nodal values of "all" variables have participated. Every wave peak, 
whether sinusoidal or sawtooth type, has been recognized. Shock waves interacting with 
turbulence, shear layers, or shedding of vortices, or effects of chemical reactions can be 
reflected in calculations of (17 a, b, c) and could have eventually contributed to the 
solution of (16) for determination of stability or instability. 

Solutions to the nonlinear ordinary differential equation (16) may be obtained using 
the fourth-order Runge-Kutta method. Iterations will continue until convergence. 

7.5 SOLUTION PROCEDURE 

It is clear that the solution consists of three parts. First the Navier— Stokes 
equations are solved. Then the results of the Navier— Stokes solutions are used to calculate 
the energy growth rate parameters. Finally the energy growth factor is computed by 
solving the nonlinear, nonisentropic stability equation. Step— by— step solution procedures 
are described as follows: 

(1) With appropriate boundary and initial conditions, solve the Navier— Stokes 

equations. Initially, the mean pressure, p, and temperature, T, based on the ideal 



28 


gas law, are specified everywhere. At the inlet for liquid propellants, however, the 
oscillatory pressure is specified [p = p + d sin art)], where d is the % disturbance 
and w is the frequency. For the solid propellants, the burning surface gas normal 
velocity (v) is computed from the response function and propellant burning rate 
such that v = v(l 4- d sin uA). 

(2) Calculate p, v^ p, and T. Taylor— Galerkin finite element method with adaptive 
meshes is used in Navier— Stokes solutions. 

(3) Advance time steps (At) of Navier— Stokes solutions to obtain wave oscillations to 
cover at least one wave period. At is determined continuously which satisfies an 
acceptable Courant number. 

(4) Take time averages for the period nAt with n chosen such that at least one peak 
wave period is covered. These time averages lead to p, Vj, and p. 

(5) Calculate the fluctuation quantities as p' = p — p, v- = — v i} p' = p — p, where p, 

v i} and p represent Navier— Stokes solutions. 

(6) Calculate the energy growth rate parameters a { , a 2 , and a 3 from (17a, b, c) using 
the results of step 5, above. 

(7) Solve the nonlinear differential equation (16) using the Runge— Kutta method with 
the initial condition e = 1 at t = 0, corresponding to neutral stability. 

(8) Repeat steps 1 through 7 until the desired length of time has been advanced. 

Note that for each time— average period in step 4, above, instability and stability are 

determined by e > 1 and e < 1, respectively, with e = 1 being the neutral stability. If the 

system is found to be unstable, it is not necessary to proceed to the next time step. 

However, for the entire ranges of time for which Navier— Stokes solutions are available, the 

stability analysis may be performed, even if instability has been found in previous time 



29 


steps. This is so because Navier— Stokes solutions are independent of the stability analysis 
as formulated here. Rather, the stability analysis here determines the state of stability or 
instability based on current flowfield as calculated from the Navier-Stokes solution. The 
following applications are based on computer code ECI— 2. 

7.6. APPLICATIONS 

7.6.1 LIQUID PROPELLANTS 

Case 1 Laminar Compressible Nonreacting Flows 
Figure 2a shows the axisymmetric geometry of combustion/thrust chamber with 965 
triangular elements and 528 nodes as a consequence of adaptive mesh process. Notice that 
in the vicinity of the walls approximately 80% of the nodes are concentrated resulting in 
prominent boundary layers. To demonstrate the capability of the code, initially, the 
steady state flow field without disturbances is examined. The velocity vectors and 
contours of Mach number, pressure, and temperature are shown in Fig. 2b, c, d, and e, 
respectively. Formation of boundary layers (Fig. 2b), separation of boundary layers and 
weak shock waves downstream (Fig. 2c, d), and decrease of temperature downstream and 
toward the wall (Fig. 2e) are evident. 

With disturbances of d = 10%, 20% and 30% imposed on the mean pressures at the 
inlet, the Navier-Stokes transient analyses are performed. The time steps are continuously 
adjusted to satisfy acceptable Courant numbers, 0.2 ^ CN ( 0.4, for convergence. The 
graphical representation of oscillations of all variables at every node versus time is 
overwhelmingly complex. Therefore, wave forms only for d = 30% and p = 3000 psi at 
three selected positions, A at (1.8, 11.43 cm), B at (31.75, 11.43 cm), and C at (63.5, 6.55 
cm) are shown in Fig. 3. Note that, although the sinusoidal input is provided at the inlet, 
the oscillations downstream become nonlinear, possibly of sawtooth type. 
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In Fig. 4, the energy growth factors £ for p = 500 psi are shown for various % 
disturbances. It is seen that for d = 10%, the energy growth factor remains in the stable 
region £ < 1. The dotted lines and solid lines indicate the results for the linear (Eq. 20) 
and nonlinear (16) cases, respectively. As the disturbance increases (d = 20%) the energy 
growth factor reaches the neutral stability, e = 1 at t ~ 0.015 sec. For d = 30%, the energy 
growth factor increases further (e = 1.038) at t * 0.2 sec. Notice that the linear analysis 
underestimates the stability if stable whereas it underestimates the instability if unstable. 
The general trend is that instability is proportional to the percent disturbances. 

As the mean pressure increases (p = 3000 psi), the possibility of instability increases 
as shown in Fig. 5, with the conclusion that instability is proportional to the mean 
pressure. Recall that for d = 10%, p = 500 psi, stability prevailed throughout whereas 
with d = 10%, p = 3000 psi neutral stability has been reached. The peak values of e for p 
= 3000 psi are significantly larger than those for p = 500 psi. 

Case 2 Turbulen t Compressible Nonreacting Flows 
The discretized geometry for a steady state turbulent compressible nonreacting flow 
is shown in Fig. 6a with a total of 2688 elements and 1416 nodes. As seen in Fig. 6b 
(velocity vectors) and Fig. 6c (Mach number contours), the boundary layers are thinner 
than in the case of laminar flow, leading to turbulent shock wave interactions. It is clear 

that gradients of pressure (Fig. 6d) and temperature (Fig. 6e) are larger than in the case of 
laminar flow. 

Wave forms of transient turbulent nonreacting flow with disturbances, d = 30%, 
p = 3000 psi, at the three positions are shown in Fig. 7. Although the wave forms in 
turbulence at these positions do not seem much different from the case of laminar flow, it is 
quite possible that wave forms at other locations where turbulent velocity gradients are 
significant would be drastically changed in contrast to the laminar flow. 
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The energy growth factors for turbulent flow are shown in Fig. 8 for p = 500 psi and 
Fig. 9 for p = 3000 psi. It is interesting to see that the effect of turbulence is to increase 
instability as compared with the laminar flow. The general trend other than the turbulent 
effect, however, remains the same as the laminar flow. That is, instability is proportional 
to disturbances and the mean pressure. The Unear analysis again shows that stabiUty is 
underestimated if stable and instabiUty is underestimated if unstable. 

Case 3 rhPTnirallv ^parting Laminar Compressible Flow s with Hy drogen/Oxyge n 

Combustion 

Figure 10 shows the discretized geometry for a steady state chemically reacting flow 
without disturbances. A total of 1580 elements and 844 nodes are used. The chemical 
reactions considered are shown in Appendix B. The mass fractions at the inlet are 0.111 
for H 2 and 0.889 for 0 2 , and the inlet velocity is 500 m/s. In this analysis the effect of 
viscosity is ignored to ensure an enhanced computational convergence, which leads to a 
flow without boundary layers along the wall. Due to the finite rate chemistry and stiffness 
arising from the chemical source terms, the computational convergence is rather slow. In 
the region of chemical reactions the velocity is decreased (Fig. 10b) and the Mach number 
contours are spaced widely apart (Fig. 10c), resulting in a rapid increase of pressure (Fig. 
lOd), but no evidence of shock discontinuities (Fig. 10c, d) between the throat and the 
downstream nozzle area. Temperature increases in the region of combustion in a sharp 

contrast to the nonreacting cases (Fig. 2e and Fig. 6e). 

The contours of 8 species are shown in Fig. 11, beginning with the reactants, 
hydrogen (Fig. 11a) and oxygen (Fig. lib), followed by products, H (Fig. 11c), H0 2 (Fig. 
lid), H 2 0 (Fig. lie), H 2 0 2 (Fig. Ilf), O (Fig. llg), and OH (Fig. llh). It is observed that 
H 2 is depleted halfway between the inlet and throat whereas 0 2 prevails further 
downstream before it is depleted behind the throat. Most of the products (including the 
radicals), H, H0 2 , O, and OH rapidly increase from the inlet and become maximum as 
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they pass through the throat. In contrast, H 2 0 and H 2 0 2 gradually increase downstream 
with a constant rate. 

Figure 12 shows the wave forms of the transient chemically reacting flow with 
disturbances, d = 30% and p = 500 psi, at the three positions considered earlier. Notice 
that, with chemical reactions, the frequencies of waves are very low and the response is 
quite slow particularly at downstream locations. Once again, the results plotted in Fig. 12 
are misleading because oscillations in other locations (over 800 nodes) may prove to be 
significantly different. After all, combustion instability is determined by the stability 
equation, not by the appearance of oscillations observed at random locations. 

In Fig. 13, the energy growth factors for the chemically reacting flow are examined. 
First of all, the linear analysis (dotted line) shows an apparent faulty prediction of 
instability for d = 10%, in which the nonlinear analysis indicates e ~ 0 throughout the time 
segment investigated. Obviously, this is an indication of unreliability of the linear 
analysis. For d = 20%, however, the linear analysis appears to give the consistent results 
similar to the earlier examples in that the linear analysis underestimates the stability when 
stable and underestimates the instability when unstable. The nonlinear analysis shows 
that, for d = 20% and d = 30%, the peak values of energy growth factors are significantly 
larger than the nonreacting cases. Does this imply that chemically reacting flows always 
tend toward instability? An affirmative answer to this question is premature. In fact the 
entire investigation presented here is subject to the future verification by experimental 
measurements. 

7.6.2 SOLID PROPELLANTS 

As an example we select the motor data used at Naval Weapons Laboratory (NWC, 
Motor #9 in [152]). The following propellant and gas data are used: propellant density p s 
= 1689 kg/m 3 ; propellant burning rate v 0 = 0.24 in/s; response function R = 0.85 
(frequency u = 300 Hz) as determined from the NWC test data (see Fig. 14); burning 
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surface pressure p - 1000 psi; gas density p = 7.1 kg/m 3 ; burning surface normal velocity v 

= 0.784 m/s; gas temperature T = 2883' R; Reynolds number Re = 10,000; viscosity p = 
1.4138 x 10' 5 kg/m/s. 

For calculations required for the stability analysis, we use a smaller motor (L/D = 
6) in order to reduce computer time with the geometry as shown in Fig. 15. The finite 
element intermediate and final adaptive meshes are shown in Fig. 16a and Fig. 16b, 
respectively. Pressure and velocity data at the burning surface are shown in Fig. 17a and 
Fig. 17b, respectively. 

Fig. 18 shows the streamline contours 0.0 < ip < 980 with Aip = 2.8 as they reach the 
steady state condition. Here the Navier-Stokes time step is chosen to satisfy the Courant 

number equal to approximately 0.4. Steep velocity gradients are apparent in the vicinity 
of the nozzle. 

Figs. 19 through 26 indicate wave motions at locations A and B (Fig. 15) for 
pressure, longitudinal velocity, and radial velocity. Note that, as disturbances increase 
from 10% to 80% at A, the frequency and amplitude for the pressure increase significantly 
as seen in Figs. 19a through 19d. The longitudinal velocity turns to steep-fronted 
N-waves as disturbances increase (Figs. 20c through 20d). Similar trends occur in the 
radial velocity as shown in Figs. 22a through 22d. At location B, the pressure peak is 
smaller than at A (Figs. 22a through 22d). The longitudinal velocity increases drastically 
and there is a clear evidence of shock waves as demonstrated by the steep-fronted wave 
forms as disturbances increase (Figs. 23a through 23d). The radial velocity, however, 
decreases significantly with waveforms being quite irregular particularly at large 
disturbances (Figs. 24a through 24d). 

Figures 25 and 26 represent waterfall plots for pressure and longitudinal velocity 
modes, extended to time t = 0.12 sec. Initial pressure disturbances upon ignition disappear 
as time progresses, but as disturbances increase, wave motions become prominent at time 
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t = 0.05 sec. with the frequency of approximately 990 Hz (d = 50%, 80%) as shown in Figs. 
25a through d. For the longitudinal velocity (Figs. 26a through d), however, a prominent 
peaks occur at d = 30% at approximately 980 Hz but are gradually reduced as disturbances 
increase to 80%, but instead smaller peaks appear at 1300 Hz, 2500 Hz, 2500 Hz, and 3100 
Hz. These phenomena correspond to the stability calculations as discussed below. 

Finally, the energy growth factors for various % disturbances are shown in Figs. 27a 
through 27d. Solid and dotted lines represent the nonlinear (Eq. (16)) and linear (Eq. 
(20)) analyses, respectively. It is seen that stability condition (0 < c < 1) is maintained for 
20% < d < 50% but instability (e > 1) occurs at t = 0.03 sec. for d = 80%. This 
corresponds to the waterfall peak at 990 Hz which begins at t = 0.03 sec. Recall, however, 
that the energy growth factor which determines instability represents the entire system 
behavior contributed from oscillations of all variables. It is interesting to note that the 
linear analysis always underestimates the stability if stable but underestimates the 
instability if unstable. It should be remarked that, for this short motor (L/D = 6), it took 
the large % disturbance (d = 80%) to cause the motor unstable. In a separate analysis 
with the full length motor (NWC #9 motor, L/D = 35.6) the motor became unstable at d 
= 20%. 


7.7 CONCLUSIONS 

The entropy controlled instability method has been applied to various problems in 
laminar flows, turbulent flows, and reacting flows for determination of stability conditions. 
The following conclusions are reached: 

(1) Instability increases with an increase of disturbances. 

(2) Instability increases with an increase of the mean pressure. 

(3) Instability increases as laminar flows are changed to turbulent flows. 

(4) Instability due to production of radicals under the finite rate chemistry is significant 


for the case investigated. 
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(5) The linear stability analysis underestimates stability if stable and underestimates 
instability if unstable. 

(6) The correct stability analysis calls for at least the third order nonlinearity. 

(7) Navier— Stokes solutions are capable of simulating the previously observed physical 
phenomena such as pressure and velocity coupling, flow turning, DC shifts, distributed 
combustion, pulsing, limit cycles, triggering, etc. 

(8) Effects of all physical phenomena are then reflected in the energy growth rate 
parameters 04, a 2 , and a 3 and subsequently in the energy growth factor e. 

(9) Future studies are required to validate the present theory in comparison with 
experimental results for those cases other than examined in this paper. 

8. CONCLUDING REMARKS 

Unstructured grids, adaptive methods, and vector processing have been applied in 
developing a major computer code for the solution of reacting flow Navier— Stokes system 
and determination of combustion instabilities. Implementation of parallel processing is in 
progress and preliminary results are expected shortly. 

This report does not include evaluation of codes developed by other organizations. 
Instead, the basic criteria for accuracy and efficiency have been established, and some 
applications have been made on rocket combustion problems. Research toward an ultimate 
goal, the most accurate and efficient CFD code, is in progress and will continue for years to 


come. 
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APPENDIX B 
REACTION EQUATIONS 


H 2 + 0 2 — 20H 

h + o 2 = oh + o 

OH + H 2 = H 2 0 + H 
O + H 2 = OH + H 
20H = H 2 0 + 0 
H + OH + M = H 2 0 + M 
2H + M = H 2 + M 
H + 0 2 + M = H0 2 + M 
H0 2 + OH = H 2 0 + 0 2 


H0 2 + H — H 2 + 0 2 
H0 2 + H = 20H 
H0 2 + O = OH + 0 2 
2H0 2 = H 2 0 2 + 0 2 
H0 2 + H 2 = H 2 0 2 + H 
H 2 0 2 + OH = H0 2 + H 2 0 

h 2 o 2 + H = oh + ho 2 
h 2 o 2 + o = oh + ho 2 

H 2 0 2 + M = 20H + M 



Any variable (f) 
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(a) Concept of ECI method 
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Fig. I Schematic overview of ECI method 
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Fig. 3 


l^ar nrareacting flow with disturbances. Procure 

** timc ■* ▼anons positions. Position 
A (1.8, 11.43 cm), position B (31.75, 11.45 cm), position C (63.5, 6.55 cm) 


52 


a 

a 

a 



a 

a 



QD 

ro 



TIME [ s©c ) 


Line** Nonlinear 


Fie. 4 Energy growth factor vs time, transient laminar flow with disturbances d = sno 
psi, f (inlet) = 1022 Hz, M (inlet) = 0.2 ’ P 
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Fig. 8 Energy growth factor vs time, transient turbulent flow with disturbances, 
p = 500 psi, f (inlet) = 1022 Hz, M (inlet) = 0.2 
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Fig. 9 Energy growth factor vs time, transient turbulent flow with disturbances, 
p = 3000 psi, f (inlet) = 1022 Hz, M (inlet) = 0.2 
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Frequency — Hz 


Fig 14 NWC response function measurement, Reference 15 , used for determination of 
burning surface normal velocity 
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Fig 15 Motor geometry for stability calculations, A (x = 5.49*, y = 0.79’), 
B (x = 13’, y = 0.84’) 



Fig 16a Finite element intermediate adaptive mesb confi miration 
332 elements, 219 nodes ’ 



Fig 16b 


Finite element final adaptive mesh configuration, 4040 elements, 2159 nodes 
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Burning surface mean pressure p, oscillatory motion mo d e l e d by 
p’ = p (1 + dsin ut)j d = percent disturbance, u = driving frequency 
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Burning surface mean vdodty v, oscillatory motion calculated 
from response function 


Fig 18 Streamline contours, minimum ^ = 0, mtmmnm $ = ggo with = 2.8 
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Fig 19 Pressure waveforms at location A for various % disturbances 
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Fig 20 Longitudinal velocity waveform, at location A for various % disturbances 



^8 21 Radial velocity waveform)! at location A for various % disturbances 
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Pressure waveforms at location B for various % disturbances 
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Rg 23 L0 " gUudl “ 1 ' vclodt » waveforms at location B for various % disturbances 
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(a) d = 10% 
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Frequency (Hz) 
(c) d = 50% 



(b) d = 30% (d) d = 80% 


Waterfall plots for pressure at location A for various % disturbances 
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Fig 25 



(b) d = 30% 


(d) d = 80% 


Fig 26 


Waterfall plots for longitudinal velocity at location A for 
various % disturbances 
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